function D = view_spatial_confounds(S);

% D = view_spatial_confounds(S);
%
% S.D        = spm eeg dataset containing entries in sconfounds (required)
%  .method   = 'view' (just show confound topographies - default) or 
%              'reject' (GUI-based confound rejection)
%
% TO PLOT COMPONENT TIMECOURSES ALSO INCLUDE:
% S.Draw     = 'raw' dataset from which spatial confounds will be removed
%              (can be the same as S.D)
%  .sampwin  = timewindow (in samples) of raw dataset to plot
%  .compchan = index of 'artifact channel' (EOG/ECG) to plot alongside
%              component timecourses

% LH 030809

f1 = sfigure;

D = S.D;
try
    Draw = S.Draw; %raw data structure, needed to plot timecourses
    sampwin = S.sampwin; %sample window, needed to plot timecourses
    hasraw = 1;
    rawdat = Draw(Draw.meegchannels,sampwin,:);
    f2 = sfigure;
    if isfield(S,'compchan')
        sfigure;
        ccdat = Draw(S.compchan,sampwin,1);
        plot(ccdat);
        title('Artifact channel (for comparison)');
        hascc=1;
    else
        hascc=0;
    end
catch
    hasraw = 0;
    warning('Problem accessing raw data; component timecourses will not be shown.')
end

try
    method = lower(S.method);
catch
    method = 'view';
end


if any(strncmp(D.chantype,'MEG',3))
    hasmeg=1;
    sens    = D.sensors('MEG');
    laycfg = [];
    laycfg.grad = sens;
    laymeg = ft_prepare_layout(laycfg);

else
    hasmeg=0;
end
if any(strncmp(D.chantype,'EEG',3))
    haseeg=1;
    sens    = D.sensors('EEG');
    laycfg = [];
    laycfg.elec = sens;
    laycfg.rotate = 0;
    layeeg = ft_prepare_layout(laycfg);
else
    haseeg=0;
end
comp = D.sconfounds;
for i = 1:size(comp,2)
    sfigure(f1);clf;
    if hasmeg
        subplot(221);topoplot([],laymeg.pos(1:102,1),laymeg.pos(1:102,2),comp(1:3:306,i),laymeg.label(1:102));
        subplot(222);topoplot([],laymeg.pos(1:102,1),laymeg.pos(1:102,2),comp(2:3:306,i),laymeg.label(1:102));
        subplot(223);topoplot([],laymeg.pos(1:102,1),laymeg.pos(1:102,2),comp(3:3:306,i),laymeg.label(1:102));
    end
    if haseeg
        if hasmeg;subplot(224);end
        topoplot([],layeeg.pos(1:60,1),layeeg.pos(1:60,2),comp(307:366,i),layeeg.label(1:60));
    end
    if hasraw
        sfigure(f2);clf;
        confound = comp(:,i);
        plot(pinv(confound)*rawdat);
        title('Timecourse of component');
        if hascc
            correlation = corrcoef(ccdat,pinv(confound)*rawdat)
        end
    end
    if strcmp(method,'reject')
        sfigure(f1);
        if(S.auto_remove>0)
            if(i<=S.auto_remove)
                keepme{i}=1;
                [S.fname '_eb_spatial_comp' num2str(i)]
                print('-dpng',[S.fname '_spatial_comp' num2str(i)]);
            end;
        else
            keepme{i} = spm_input('Keep this component as confound?','+1','Yes|No', strvcat('Y', 'N'));
        end;
    end
end

if strcmp(method,'reject')
    rejectme = strcmp(keepme,'N');
    comp(:,find(rejectme))=[];
    ncomp = size(comp,2);

    [sel1, sel2] = spm_match_str(D.chanlabels(D.meegchannels), D.chanlabels(setdiff(D.meegchannels, D.badchannels)));
    sconf = [];
    sconf.label = D.chanlabels(D.meegchannels);
    sconf.coeff = nan(length(sconf.label), ncomp);
    sconf.coeff(sel1, :) = comp(sel2, 1:ncomp);
    sconf.bad = ones(length(sconf.label), 1);
    sconf.bad(sel1, :) = 0;
    D = sconfounds(D, sconf);
end
